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ABSTRACT 

Iterative Filtering (IF) is an alternative technique to the Em¬ 
pirical Mode Decomposition (EMD) algorithm for the de¬ 
composition of non-stationary and non-linear signals. Re¬ 
cently in [2] IF has been proved to be convergent for any 
signal and its stability has been also showed through exam¬ 
ples. Furthermore in [2] the so called Fokker-Planck (FP) fil¬ 
ters have been introduced. They are smooth at every point and 
have compact supports. Based on those results, in this paper 
we introduce the Multidimensional Iterative Filtering (MIF) 
technique for the decomposition and time-frequency analysis 
of non-stationary high-dimensional signals. And we present 
the extension of FP filters to higher dimensions. We illustrate 
the promising performance of MIF algorithm, equipped with 
high-dimensional FP filters, when applied to the decomposi¬ 
tion of 2D signals. 

Index Terms — Iterative Filtering, Multidimensional It¬ 
erative Filtering, Empirical Mode Decomposition, non-linear 
and non-stationary multidimensional signals, Fokker-Planck 
filters 

1. INTRODUCTION 

Given a non-stationary signal a hard problem is how to per¬ 
form a time-frequency analysis in order to unravel its hidden 
features. The problem becomes even harder if we want to 
handle a signal that have dimension higher than one. Such 
kind of problems are ubiquitous in real life and their solutions 
can help shedding lights in many research fields like, for in¬ 
stance, in the non-destructive detection of structural damages 
in buildings and machineries, in the quantitative canvas weave 
analysis in the art investigations of paintings, in the identi¬ 
fication of harmful airborne chemical particles by means of 
hyperspectral images analysis, in the image enhancement for 
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medical applications, in the extraction of information from 
atomic crystal images, et cetera. 

Time-frequency analysis have been substantially studied 
in the past [4] and, traditionally, linear techniques like Fourier 
spectral analysis or wavelet transforms have been commonly 
used to decompose signals, even non-stationary ones, into 
simple and stationary components and then perform on them 
a frequency analysis. These approaches, even if easy to im¬ 
plement, have limitations. In particular, they work well when 
the given signal is periodic and stationary, whereas they can¬ 
not handle properly non-stationary signals. Furthermore all 
these techniques use predetermined bases and they are not de¬ 
signed to be data-adaptive. Hence there is the need to develop 
new methods designed to handle specifically non-linear and 
non-stationary signals. 

The research on these new techniques started in 1998 
when Huang and his research group developed and released 
the very first method of this kind, the so called Empirical 
Mode Decomposition (EMD) algorithm [9]. The idea behind 
this method is to generalize the approach of classical linear 
techniques: decompose a signal into simple components and 
then perform a time-frequency analysis on each component 
separately. For this reason Huang and his group defined the 
so called Intrinsic Mode Functions (IMFs) which are simple 
functions with the following properties: their number of ex¬ 
trema is either equal to the number of zero crossings or they 
differ at most by one, and at any point their moving average 
is zero. Furthermore they developed an iterative technique, 
called the sifting process, to decompose a signal into such 
Intrinsic Mode Functions with the final goal of computing the 
instantaneous frequency [1,4] of these simple components. 

The sifting process of the EMD is structured in the fol¬ 
lowing way. Let L be an operator capturing the moving av¬ 
erage of a signal /, and let S{f){t) = f{t) — be its 

fiuctuation part. Then the first IMF is given by IMFi(f) = 
lim^^oo The limit is such that the moving av¬ 

erage of IMFi(t) is the zero function. Assuming k > 1 
IMFs have been computed, then the subsequent IMFs are pro¬ 
duced applying the aforementioned procedure to the residual 
r{t) = f{t) — X]j=i IMFj(t). The method stops when the 



residual r becomes a trend signal. So in the end the given 
dataset is decomposed as f{t) = 

this algorithm the moving average L{f){t) is given by the 
mean function of an upper and a lower envelope, where upper 
and lower envelopes are cubic splines connecting local max¬ 
ima and local minima of f{t), respectively. This method has 
received a lot of attention in the last decade. In fact, many 
researchers have applied it to solve several open problems in 
diverse research fields. 

We point out that there are known examples showing that 
this technique can be unstable to small perturbations of the 
initial data. To overcome this issue Huang et al. developed 
the Ensemble Empirical Mode Decomposition (EEMD) [17] 
where the IMEs are taken as the mean of many different trials 
produced using the EMD algorithm. Since cubic splines are 
used repeatedly in the iterations of both EEMD and EMD, 
a rigorous proof of the convergence of these methods is still 
missing. 

Inspired by the seminal work by Huang et al., in the 
last years many research groups started working to the de¬ 
velopment of alternative techniques to the EMD algorithm. 
Eor instance Daubechies and her research group devised 
the Synchrosqueezed wavelet transform [5], Gilles created 
the Empirical wavelet transform [7], Hou et al., using the 
multicomponent amplitude modulation and frequency mod¬ 
ulation (AM-EM) representation [15], developed the Sparse 
time-frequency representation method [8], and many oth¬ 
ers [6,12, 13,16]. All these alternative methods make use 
of optimization techniques for the decomposition of a given 
non-stationary and non-linear signal. These approaches 
require the a priori selection of a suitable basis for the decom¬ 
position. 

In 2009 Zhou et al. developed another alternative algo¬ 
rithm called Iterative Eiltering [10]. This method, as the orig¬ 
inal EMD, has an iterative structure and, unlike all the tech¬ 
niques based on optimization, it does not require any initial 
assumption on the signal. Therefore it is able to produce de¬ 
compositions that are completely data driven. This method 
has the same structure of the sifting process, but now the mov¬ 
ing average L{f){t) is derived by convolution of the given 
signal f{t) with filters like, for instance, a double average 
filter. This new way of computing moving averages allows 
to perform a rigorous analysis of this technique. In particu¬ 
lar, under mild sufficient conditions on the filters used in the 
convolutions, the analytic convergence of Iterative Eiltering 
applied to generic signals is ensured and the components 
produced in the decomposition of a signal can be derived from 
an explicit analytic formula [2]. 

Some of the aforementioned methods, initially developed 
to decompose ID datasets, have been already generalized to 
handle two and higher dimensional signals, like for instance 
the multidimensional ensemble empirical mode decomposi¬ 
tion method [18] or the Synchrosqueezed wave packet and 
curvelet transform [19]. However this is not the case of Iter- 


MIF Algorithm IMFs = MIF(/) 

IMFs = {} 

while the average number of extrema of / > 2 do 
compute the filter support Q for / 

/i = / 

while the stopping criterion is not satisfied do 
/„+i(x) = /„(x) - f„(x + t)w(t)d''t 
n = n -h 1 

end while 

IMFs = IMFs U{/n} 

f = f - fn 

end while 

IMFs = IMFs U {/} 


Table 1. Multidimensional Iterative Filtering pseudocode 

ative Filtering. Therefore in this paper we address the prob¬ 
lem of extending Iterative Filtering and the so called Fokker- 
Planck (FP) filters [2] to higher dimensions. We introduce, in 
particular, the Multidimensional Iterative Filtering algorithm 
and the Generalized Fokker-Planck (GFP) filters. We show 
also the performance of this method, equipped with GFP fil¬ 
ters, on artificial and real life signals. 

2. MULTIDIMENSIONAL ITERATIVE FILTERING 
ALGORITHM AND GENERALIZED 
FOKKER-PLANCK FILTERS 

The Multidimensional Iterative Filtering (MIF) technique 
consists of an Inner and an Outer Loop. In the inner loop the 
method computes an IMF of a given /c-dimensional signal 
/ as the limit of the sequence generated by subtracting from 
the signal its moving average iteratively. In the outer loop we 
simply update the signal by removing from it the previously 
computed IMFs. The outer loop is iterated until the remain¬ 
der becomes a trend signal. Whereas the inner loop should 
be theoretically iterated until the moving average becomes 
a zero function. However, in the numerical implementation 
of the inner loop, we use a stopping criterion to discontinue 
the iterations. The pseudocode of this algorithm is given in 
Table 1, where / is the /c-dimensional signal we want to 
decompose and re G R* is a filter function with finite support 
n c R^. 

As a filter w we need to have, first of all, nonnegative 
functions which are finitely supported on Q C with 
w(t)d^t = 1. However not every w function which ful¬ 
fils the previous properties is going to ensure the convergence 
of this technique. A well chosen class of filters has to be 
adopted. From what is known for the ID Iterative Filtering 
algorithm [2], we can conjecture that in higher dimensions a 
class of filters which guarantees the MIF convergence is given 



by a proper extension of the ID Fokker-Planck (FP) filters. 
Such ID filters, derived from the solution of Fokker-Planck 
partial differential equations [2], have the nice property of be¬ 
ing extremely smooth since they are infinitely differentiable 
at any point of their domain. To produce axial symmetric 
Fokker-Planck filters of dimension higher than one there is 
no need to numerically solve high-dimensional partial differ¬ 
ential equations. We can make use of highly accurate numer¬ 
ical solutions in ID to generate, first, numerical 2D Fokker- 
Planck filters and then higher dimensional ones. This can 
be achieved by simply scaling, resampling, and then rotating 
around one axis a ID filter many times to produce a 2D ver¬ 
sion of it. Higher dimensional filters can be produced simply 
iterating this procedure. We call them Generalized Fokker- 
Planck (GFP) filters. We observe that this approach allows to 
produce also GFP filters with an ellipsoidal support D. 

Once a filter shape has been selected it can be used for 
every dataset. The problem is, in order to extract meaningful 
IMFs, how to select a proper support Ft for the filter simply 
based on the signal we want to decompose. Following what 
proposed in [10] for the ID case, we compute for each dimen¬ 
sion of the signal the average ID length of the support. Then 
we use this information to either find the radius of a spherical 
support or to identify the radii of an ellipsoidal Ft C . 

About the stopping criterion, there are many possibile op¬ 
tions. One way of doing it, as mentioned in [2,9] for the ID 
case, is to consider the relative change in and discontinue 
the inner loop as soon as a prefixed threshold is reached. 

Finally, we observe that this extension to higher dimen¬ 
sions can be applied in the same way also to the so called 
Adaptive Local Iterative Filtering (ALIF) method [2], which 
is a ID generalization of the Iterative Filtering algorithm. 


3. NUMERICAL EXAMPLES 

In the following we show the performance of the Multidi¬ 
mensional Iterative Filtering (MIF) algorithm, equipped with 
Generalized Fokker-Planck (GFP) filters, applied to both ar¬ 
tificial and real life signals. For simplicity from now on we 
consider only bidimensional signals, however we point out 
that the MIF method can handle datasets of any dimension. 


3.1. Example 1 

We start considering the artificial signal showed in Figure 1 
containing the mixture of two non-stationary sinusoidal 
damped signals. 

If we apply the MIF algorithm to this signal we obtain the 
two IMFs plotted in Figure 2. 

To have a better understanding of the accuracy of the de¬ 
composition we can plot sections of these IMFs taken along 
anti-diagonals. This approach allows to easily compare the 
IMFs with the ground truth, as showed in Figure 3. 
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Fig. 1. Artificial bidimensional dataset containing a mixture 
of two non-stationary sinusoidal signals. 
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Fig. 2. IMFs produced using MIF. 
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Fig. 3. Sections along the anti-diagonal of the IMFs. 



Fig. 4. A non-stationary and non-smooth signal and its mid¬ 
dle vertical section. 
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Fig. 5. First (left) and second (right) IMF produced by MIF 



Fig. 6. Middle vertical 
spending ground truth. 
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3.2. Example 2 

For the second example we consider the signal of Figure 4 
which is the mixture of a smooth sinusoidal signal with a non¬ 
smooth and non-stationary one. If we run the method, using 
a 2D generalized Fokker-Planck filter, we can produce a first 
IMF which is a non-smooth and non-stationary signal and 
a second one which is smooth and stationary, as showed in 
Figure 5. To have a better understanding of the kind of sig¬ 
nal we are dealing with and to appreciate the performance of 
the decomposition algorithm, we plot middle vertical slices 
of the original signal, IMFs, remainder and ground truth in 
Figures 4, 6 and 7. This example shows the ability of this 
method to separate components completely different in na¬ 
ture, like smooth and non-smooth IMFs, using a fixed filter 
function. It would be hard to produce similar results with any 
optimization technique. 



Fig. 7. The Remainder, its middle vertical section and the 
corresponding ground truth. 



Fig. 8. Contrast-enhanced spectral-mean image of the hyper¬ 
cube dataset. 



Fig. 9. Pixel classification using the raw (left) and pre- and 
post-processed (right) dataset. 


3.3. Example 3 

As a real life signal we consider the case of a hyperspectral 
image F G x v x ^ frequency channels, h x v pix¬ 

els pij and their corresponding signatures Sij G R^. In Fig¬ 
ure 8 we show the hyperspectral image of an area where some 
chemical has been released in the atmosphere. The problem 
is, given the spectral signature of a chemical Sc G R^, to clas¬ 
sify the pixels of the hypercube in order to identify the chemi¬ 
cal position in the air [11]. This can be done using a classifier 
like the Adaptive Cosine Estimator (ACE) which assign to 
each pixel a value as follows: 


y{Pij)ACE 






( 1 ) 


where S G R^ ^ ^ is the covariance matrix obtained consid¬ 
ering each pixel as an observation and each frequency as a 
variable. 

If we apply (1) to the raw hypercube data we obtain the 
classification showed in Figure 9 (left), where the darker the 
pixel is the higher is the probability that it contains the chem¬ 
ical, while the lighter the lower is the probability. We observe 
that many isolated pixels are classified in black, those are all 
false alarms that can be due to various factors like noise or 
sensor malfunction. 

We can use the MIF algorithm to first pre-process the hy¬ 
percube, removing from each frequency channel the first IMF, 
and then we reapply MIF to post-process the ACE classifica¬ 
tion values, removing the first IMF produced in the decom¬ 
position. In doing so the performance of the classifier are 
improved as showed in Figure 9 (right). For more details on 
this topic we refer the interested reader to [3]. 






















































4. OUTLOOK 

Inspired by the convergence and stability properties of Iter¬ 
ative Filtering and the extreme smoothness of the Fokker- 
Planck (FP) filters [2], in this paper we introduce the Multidi¬ 
mensional Iterative Filtering (MIF) algorithm for the decom¬ 
position of multidimensional non-stationary signals. Further¬ 
more we extend FP filters to higher dimensions and we test the 
MIF technique, equipped with these filters, on both artificial 
and real life datasets. These tests show the ability of MIF to 
properly decompose a non-stationary signal into IMFs even 
of completely different nature. 

This is just a preliminary work in the analysis of the MIF 
method. There are several open problems that need to be ad¬ 
dressed like, for instance: Is there a simple and effective way 
to identify and use higher dimensional position of minima and 
maxima of a multidimensional signal in the derivation of the 
filter support? What about the convergence of MIF? Are there 
sufficient conditions on the higher dimensional filters that en¬ 
sure the convergence of this method? What are the connec¬ 
tions with previously developed filtering techniques like the 
iterative steering kernel regression method [14]? 

Finally we point out that once a multidimensional signal 
is decomposed into IMFs the next step would be to apply a 
time-frequency analysis on each of them. Currently the most 
common way to perform a time-frequency analysis on ID 
IMFs is via the estimation of instantaneous frequencies by 
means of the so called Hilbert Transform [1,4,9]. However 
this approach cannot be directly extended to handle higher di¬ 
mensional IMFs, so it remains an open problem how to com¬ 
pute instantaneous frequency directly in higher dimensions. 
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